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ABSTRACT 

We describe a simple efficient algorithm that allows one to construct Monte-Carlo re- 
alizations of merger histories of dark matter halos. The algorithm is motivated by the 
excursion set model (Bond et al. 1991) for the conditional and unconditional halo mass 
functions. The forest of trees constructed using this algorithm depends on the under- 
lying power spectrum. For Poisson or white-noise initial power-spectra, the forest has 
exactly the same properties as the ensemble of trees described by Sheth (1996) and 
Sheth & Pitman (1997). In this case, many ensemble averaged higher order statistics 
of the tree distribution can be computed analytically. For Gaussian initial conditions 
with more general power-spectra, mean properties of the ensemble closely resemble the 
mean properties expected from the excursion set approach. Various statistical quanti- 
ties associated with the trees constructed using our algorithm are in good agreement 
with what is measured in numerical simulations of hierarchical gravitational clustering. 

Key words: galaxies: clustering - cosmology: theory - dark matter. 



INTRODUCTION 



It is widely believed that the massive dark matter halos 
which exist today have grown from small initially Gaussian 
density fluctuations. Many questions of astrophysical inter- 
est can be addressed using Monte-Carlo realizations of the 
merger histories of such massive dark matter halos. An effi- 
cient method for generating Monte-Carlo realizations of the 
clustering process is essential for addressing such questions 
(e.g. Kauffmann & White 1993). Ideally, such an algorithm 
should produce results that are consistent with other known 
properties of the clustering process. 

In this paper, we develop an algorithm which is mo- 
tivated by known excursion set results (Bond et al. 1991; 
Lacey & Cole 1993). For example, these earlier results pro- 
vide expressions for the mean number of halos of mass 
m that are later incorporated into a larger halo of mass 
M > m. For the special case of clustering from Poisson, or 
white-noise initial conditions, a model for the higher order 
moments of this subclump distribution exists (Sheth 1996; 
Sheth & Pitman 1997). Appendix ^ of this paper shows 
that these higher order moments are consistent with the fact 
that mutually disconnected regions in Poisson or white-noise 
Gaussian distributions are mutually independent. Section ^ 
describes an algorithm which uses this fact to generate a 
forest of merger trees. The algorithm reproduces the mean 
and higher-order Poisson and white-noise results exactly. Al- 



though the algorithm is not of binary split type, trees con- 
structed usin g bin ary splits provide good approximations to 
ours (Section 3.3). 



Unfortunately, previous work does not provide expres- 
sions for the higher order moments of the subclump distri- 
bution for other initial power spectra (e.g. Sheth & Pitman 
1997). Section ^ shows that, to describe the trees associated 
with arbitrary Gaussian initial conditions, a simple modifi- 
cation of our white-noise algorithm produces good approxi- 
mations to the known (excursion set) mean values. Section ^ 
shows that, for a range of power spectra of current interest, 
our algorithm provides higher order distributions that are 
similar to those measured in numerical simulations of hier- 
archical gravitational clustering. A final section summarizes 
our results. It argues that, in addition to allowing one to 
generate a forest of merger history trees, our results are also 
useful for studying the spatial distribution of dark matter 
halos. 



2 POISSON INITIAL CONDITIONS 

The first part of this section summarizes the relevant results 
of Sheth (1996), where details of all the arguments are given. 
The second part shows how these earlier results can be used 
to derive an algorithm for partitioning halos into subclumps. 
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The algebra for some of the arguments in this section is given 
in Appendices ^ and ^. 



2.1 Excursion set and branching process results 

The evolution of clustering from an initially Poisson distri- 
bution has been studied using the excursion set (Epstein 
1983) and cloud- in-cloud formalisms (Sheth 1995). In these 
analyses, all matter is bound into clumps, sometimes called 
halos. The probability that a clump has m associated parti- 
cles is 



77(m,6) = 



(jnhy 



(1) 



where m > 1 and < 6 < 1. Equation is known as 
the Borel distribution (Borel 1942). The probability f{m,b) 
that a randomly chosen particle belongs to an m-clump is 
given by multiplying ri{m, b) by m, and dividing by the mean 
clump size mri{m, fe) = 1/(1 — 6): 



/(m, b) = m{l — b) r]{m, b). 



(2) 



Since all matter is in clumps, and /(m, b) is the fraction of 
mass in m clumps, 



^/(m,6) = l. 



(3) 



The evolution of clustering is described by the evolution of 
the b parameter, which is the ratio of the average density 
in the universe to the critical overdensity that is required 
for an object to collapse and virialize at cosmological time 
t. That is, 



l + 4(t) 



where Sdt) oc 



1.686 

a{t) ' 



(4) 



t denotes cosmological time, and a{t) is the factor by which 
the universe has expanded since some fiducial initial epoch, 
say, ti. Notice that the critical density decreases as the uni- 
verse expands, so that 6 increases with expansion. Since 6 
increases monotonically as cosmological time increases, it 
can also be treated as a pseudo-time variable. Thus, b plays 
two roles: it can be thought of as describing a spatial den- 
sity, or it can be thought of as describing time. Therefore, 
/(m, 6) can be thought of as the fraction of the total mass 
that is in m-clumps at time b, or it can be thought of as 
the fraction of the total mass that is associated with regions 
within which the average overdensity is given by b. 

The excursion set description also allows one to com- 
pute the probability that a randomly chosen particle of a 
clump having exactly M particles at the time correspond- 
ing to &o was a member of an m-clump at the earlier epoch 
fei < bo- This probability is 



/(m,6i|M, bo) 



M 



m— 1 r 

M 



bi 



(5) 



where 1 < m < M and < bi/6o < 1 (Sheth 1995). This 
implies that the average number of progenitor m-subclumps 
identified at 6i that are later (at bo > bi) in an Af-clump is 



Although the excursion set formalism provides information 
about the mean number of progenitor subclumps of any 
given parent clump, it does not provide a complete descrip- 
tion of the entire merger history tree (cf. Fig. 1 in Sheth & 
Pitman 1997). 

Let p(n, bi\M, bo), with n = (ni, • • • ,nAf), J2^Li "■i = 
n and X/j-i j = M, denote the probability that an M- 
clump at the epoch bo was previously in n subclumps at the 
epoch bi , of which there were ni single particle subclumps, 
n2 subclumps with exactly two particles, rij subclumps with 
exactly j particles, and so on. Then the branching process 
model of Sheth (1996) yields the formula 



p{n, bi\M, bo) = 



(Afboi)"-^e-^^*ni 
viM, bo) 



n 



vU,bi)"^ 



(7) 



where boi = (bo — bi). Sheth & Pitman (1997) showed that 
this same formula follows from another construction, which 
they termed the additive coalescent. Appendix^ of this pa- 
per shows that this formula is consistent with the assump- 
tion that mutually disconnected regions within a Poisson 
distribution are mutually independent. 

Notice that equation (Q) depends only on the ratio 
bi/bo, and not on the actual values of bi or bo. This depen- 
dence on the ratio bi/bo implies that the merger histories 
of M-clumps formed at the epochs bo and b'g are related to 
each other simply by rescaling the time variable. By taking 
relevant sums over the partition probability function, one 
can show explicitly that the branching process expressions 
for the clump size distribution and the merger probabili- 
ties f(m\AI) are identical to the excursion set results (equa- 
tions |l| and 1^ . However, the partition probability function 
also allows one to compute quantities that cannot be esti- 
mated using the excursion set formalism. For example, the 
excursion set approach says that the mean number of m- 
subclumps per M-clump is (M/m) f{m\AI). The partition 
probability function allows one to compute the higher order 
moments of the rij distribution also. For example. 



nil 



,bi 



(n. - a)! (nj -/3)! 

^+'3 v''{hbi)v^U,bi) viR,B) 



M,bo 



Mboi 
Mboi 



r]{M, bo) 

-+'3-1 ,,"(i,bi)/(j,bi) 



viM, bo) 
where boi = (bo — bi) and 
RB = {M -ai- Pj)B = Mbo - aibi - /3jbi 



if Ti > 

if i? = 0, (8) 



(9) 



N{m, bi\M, bo) = (M/m) /(m, bi jM, bo) 



(6) 



When a — 1 and f3 — this is equivalent to the relation 
obtained from equation When a = /3 = 1, equation (^ 
gives an estimate of the cross-correlation between i- and j- 
subclumps that are within the same M-clump. 

2.2 The partition algorithm 

The partition algorithm described below follows from the 
fact (cf. Appendix ^ that equation (Q) is consistent with 
the assumption that mutually disconnected regions within a 
Poisson distribution are mutually independent. 

Imagine partitioning an i\f-clump up into subclumps by 
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choosing first one subclump, tlien clioosing a second from 
tfie mass tfiat remains, and so on untii aii tfie mass Ad has 
been assigned to subciumps. Start by choosing a random one 
of the M particles. Equation (|^) says that, on average, the 
chosen particle will be a member of an mi-subclump with 
probability /(mi, fei |M, bo). Since 6 is a density variable, the 
halo M can be thought of as being a region of size Vm = 
Mbo/n, where n denotes the average density. If the first 
subclump has mass mi, then it is associated with a volume 
Vmi ~ m\b\/n within Vm- The remaining R\ = {M — mi) 
particles must be distributed randomly (because the initial 
distribution is Poisson) within the remaining volume Vm — 
Vmi ■ Let 6^^^ denote the average density in this remaining 



volume, and recall (equation U) that it can be represented 
by some = 1/(1 + 5^^'). Then 



n 
bW 



M ■ 



mi 



Vm — Vmi 



and so Ri (b'^^ — bi 



Mbo 



(10) 



Now the second subclump of the (A/, feo)-halo must be cho- 
sen from among the remaining Ri — M — mi particles. 
Suppose that the probability that this second subclump 
has exactly m2 particles is given by the same rule as be- 
fore, namely, by /(m2, fei|-Ri, 6^^'). If the second subclump 
has mass m2, then there are R2 = (M — mi — m2) re- 
maining particles, they are randomly distributed within 
(Vm — Kni — Kna), where Vm^ ~ m2&i/n, and the average 
density within this remaining region is pa rametrized by some 
fe'^-*, defined analogously to equation ([L0[). Suppose that the 
third subclump has exactly ms particles with probability 
/(m3,fei|7?2,6<^'), and so on. 

Let rrii denote the number of particles in the ith sub- 
clump. Let Ri denote the mass remaining after the ith pick: 



Ri = M - ^ rrij , and _Ro = M, 



(11) 



and parametrize the average density in the volume that re- 
mains after the ith pick by Then 



Rr 



Ri 



uVm - E -=1 nVm^ ^'^^01 + -^^^1 ' 

Notice that 1/6^°' = Ro/{Mboi + Robi) = 1/bo as required, 
and that 

Mboi + R,bi 



R^{b^'^ -bi)=R, (- 



R^ 



61) = Mboi (12) 



is independent of i. In this notation, equation (^ becomes 

f{m,,bi\R,^i,b^''^^) 

- "^'^^ r,(i?._i, &(-!)) ' 

- -5 Mboi ui-i-)\ — " 1 < mi < Ri-i_ 

Rt-i r){Ri-i,by^ 

and it is 

^ ri{m,i,bi) 

77(i?,_i,6('-i)) 

The probability that after n picks all M particles have been 
assigned to a subclump (i.e., m„ = Rn-i and i?„ — 0), is 
simply the product of the individual n choices above. If m 



denotes the set (mi, ■ • ■ , m„), then 

n 

p{m\M) = Y[f{m„bi\R^-i,b^'-^^) 



(Mbo 



\n-l „-M6oi 



77(M,6o) 



where the final expression uses the fact that Ro = M, fo'^-* — 
bo, and m„ — R„-i. 

Equation (^^ describes the probability that M was par- 
titioned into these n subciumps in the order (mi, ■ • ■ ,m„). 
The probability that the sequence of picks resulted in the 
set (mi, • • • , m„), whatever the order of the rrii, is given by 
summing equation (|lj) over all n! permutations of the n 
picks rrii. Let n[m] denote this set of permutations. Then 



E 

Trfml 



n-1 „-Mbn 



p{m\M) 



(A/boi)"-^e 
^ ri{M,bo) 



since, for a given m, all terms in each permutation are the 
same, except for the product J^"_j^(mi/i?i_i). So, the sum 
over permutations becomes 



E 



p{m\M) 



{Mbo 



\n-l -Mboi 



v{M,bo) 



Y]_v(-mi,bi) 



7r[7Ti] 1 



(16) 



A little thought shows that the contribution from the terms 
in the sum is unity. Appendix p3 shows this explicitly. Thus, 



^p(mlM) 



{Mbo, 



Nn-l „-Mboi " 



V{M,bo) 



Y\_v{-mi,bi) 



(17) 



Now, if all particles are indistinguishable, then subciumps 
that have the same number of particles are indistinguish- 
able from each other, so these permutations should not be 
counted separately. We can account for this by dividing the 
product above by Jlf^i where rii denotes the number of 
i-subclumps in the partition, ^ jt.; ~ n, and ^ . i rii = M. 
Thus, the probability that Af was partitioned into n is the 
same as that given by equation (^. 

This decomposition suggests the following algorithm for 
partitioning an (A/, 6o)-halo into 61-subhalos. 

Initialize: specify bl, bO, M, and set R = M. 

Set b = M(bO-bl)/R + bl . 

Choose m with probability f (m,bl I R,b) . 

Set R = R-m. 
Iterate until R = 0. 

The list of m's so obtained gives one possible partition of 
M into 61-subclumps. Repeat to generate an ensemble of 
such partitions. The previous paragraphs show that this 
algorithm will produce partitions with the same ensemble 
properties as the branching process (equation ^). 

Notice that, other than requiring that < 61/60 < 1, 
there is no other constraint on bi/bo. This means that the 
algorithm constructs partitions with the correct ensemble 
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properties in one time step. However, one is often interested 
in following the history of a halo as it becomes partitioned 
into smaller pieces through many small time steps. It is sim- 
ple to use this partition algorithm to construct such a merger 
history tree. Given M and 60, choose h\ < ho. Use the al- 
gorithm above to partition M into &i-subclumps. Call this 
partition m. Choose ti2 <b\. Partition each subclump m of 
m into its constituent subclumps using the algorithm above, 
with the initial parameters 61 — > 62, &o ^ 61, — > and 
set R — m. Then choose 63 < 62, and so on. The results 
of Sheth & Pitman (1997) and Sheth (1998) show that this 
many-step algorithm is guaranteed to give the same ensem- 
ble of partitions as the one-step algorithm described in detail 
above. 



3 GAUSSIAN INITIAL CONDITIONS 

This section shows that the partition algorithm of the pre- 
vious section can easily be generalized to partition halos 
associated with Gaussian initial fluctuations. 

Let S{t) denote the critical density required for a spheri- 
cal perturbation to collapse by time t (cf. equation ^ . Thus, 

5 is both an overdensity variable and a pseudo-time vari- 
able, as was 6 in the previous section. Also, let a^{A{) de- 
note the variance in the initial density fluctuation field when 
smoothed over regions that contain mass M on average. The 
variance depends on M in a way that depends on the power- 
spectrum; we will only be interested in initial power-spectra 
for which S{M) decreases monotonically as M increases (e.g. 
Press & Schechter 1974). 

In the excursion set approach (Bond et al. 1991; Lacey 

6 Cole 1993), the fraction of mass that is in halos which 
have mass in the range dm about m at cosmological time t 
is 



/(m, 5) dm = ^/ ^ exp ( ) ^ 



2S 



ds; 

S ' 



(18) 



where S — a'^{m,) and 5 decreases as cosmological time in- 
creases, as described above. The average number density of 
such halos is 



n(m, S) dm — (p/m) f{m, S) dm. 



(19) 



where p denotes the average density. 

Similarly, consider a halo that has mass M at time So- 
The fraction of the mass of this halo that was in subhalos 
of mass m at the earlier time 5i > So is 

f{m,Si\M,So)dm = — , provided m<M, (20) 



3.1 White-noise initial conditions 

For white- noise, a^{M) oc 1/M. In the limit of small 5, which 
corresponds to the b 1 limit (cf. equation Stirling's ap- 
proximation for the factorials reduces all the Poisson excur- 
sion set statements (equations |l| arid B) to the corresponding 
excursion set statements (|l^ and pOp"^for n — white-noise 
initial conditions (Sheth 1995, 1996). Notice that the 61/60 
dependence of the Poisson expressions translates to a de- 
pendence on {Si — So), rather than the actual values of Si or 
So themselves. This means that f{m,Si\M,So) can be writ- 
ten as f{m\M,D), where D — {Si — So). Furthermore, to 
lowest order in S, f{in\M) can be factored similarly to (|l3|). 
That is, just as in the Poisson case, M occupies a region 
V, such that M/V — p{l + So), where p is the background 
density. Similarly, the subhalo m occupies v within V, with 
m/v — p{l + Si). Thus, the density in the remaining volume 
is 



V 



— =p(l + 5'), so {Si-S')^ 



So) 



(22) 



1 - (m/A/) 

to lowest order in the S terms. Finally, since disconnected 
volumes are mutually independent in the white-noise case, 
just as in the Poisson case, the same partition algorithm 
that worked there will also work here: 

Initialize: specify dl, dO, M, and set R = M. 

Set D = (dl-dO)/(R/M) . 

Choose m with probability f(m|R,D). 

Set R = R-m. 
Iterate until R is as small as desired. 

A simple transformation of equation ( po|) shows that to 
choose m according to f{m\R,D) one needs to generate a 
Gaussian random variable. There are efficient ways to do 
this, so this algorithm is fast. 

This algorithm allows one to generate partitions of M 
whose distributional properties are described by the white- 
noise analogue of equation (^. As before, there is no re- 
striction on the size of the 'time step': {Si — So) > can be 
set as large as desired (although the algorithm is slower for 
larger time steps, since, for large Si — So most partitions are 
into many tiny pieces). And again, constructing a merger 
history tree by embedding this partition algorithm in a loop 
over time-steps is straightforward. In either case, the algo- 
rithm generates subclump distributions whose mean values 
are exactly the same as those required by the excursion set 
approach (equation The higher order moments of the 
subclump distribution are given by applying Stirling's ap- 
proximation to the corresponding Poisson expressions, and 
taking the 5^1 limit. For example, the cross-correlation 
between mi and m2 halos that are within the same M halo 



where u = 7 s ~ c"^(m), and S = a'^{M). 

Ays o j 

Thus, the mean number of m halos within an M halo is 
N{m,Si\M,5o)dm= (^^^ f{m,5i\M,So)dm (21) 
(e.g. Lacey & Cole 1993). 



c(12|0) = N{mi,Si\M,So)N{m2,Si\M -mi,S'), (23) 
where A''(m|A/) was defined by equation ( pl|) and S' was 



defined by equation 



For white-noise c(12j0) = c{21\0). 



3.2 Arbitrary Gaussian initial conditions 

The assumption that disconnected volumes are mutually in- 
dependent is wrong when the initial conditions differ from 
white-noise. Therefore, it is not obvious that the algorithm 
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above should be used to generate partitions for arbitrary ini- 
tial conditions. On the other hand, for arbitrary initial con- 
ditions, nothing is known about the higher order moments 
of the subclump distribution; only the mean is known from 
the excursion set approach. Below, we take two different 
approaches to constructing a partition algorithm for initial 
conditions that are not white-noise. Which algorithm to use 
is determined by comparison with numerical simulations in 
section ^. 

The first approach is motivated by the observation that, 
when expressed as functions of the variance rather than the 
mass, all excursion set quantities are independent of power 
spectrum. This suggests that, with a suitable transforma- 
tion, we should be able to use the algorithm above, since it is 
known to be exact for the white-noise spectrum. We do this 
as follows. Run the algorithm assuming a white-noise power- 
spectrum. However, treat each chosen m not as a subclump 
having mass m, but as a region of volume v/V = m/A4 
containing mass m, that is populated by some number u of 
objects that all have the same mass /i. Thus, u — m//i. The 
value of fj, is got by requiring that cy^ifJ-) equals the white- 
noise value, namely o-^{i-l) = 1/m. To illustrate, suppose that 
cr2(/i) = Then V is ma'^i^i)^^'^ = m*""^'/". For white 

noise, a — 1, so 1/ = 1; thus, the region v contains exactly 
one halo as in the previous subsection. For general a, v is 
neither unity nor even integer. Nevertheless, this approach 
generates partitions which are guaranteed to have the excur- 
sion set value in the mean, but which, as a result of the fact 
that \x halos always come in groups of i/, may have unnat- 
urally high values for the higher order moments. Since this 
approach simply uses the white-noise algorithm, the higher 
order moments can be calculated analytically. For example, 
if cr^(/i) = then the cross-correlation between halos is 



c(12|0) 



(24) 



where c(12|0) is given by equation (^3|) but with the white- 
noise power spectrum. 

The second approach is to assume, with no real justi- 
fication, that the algorithm of the previous subsection can 
be used directly. In this case, the ensemble of partitions will 
differ for different power-spectra only because (t^(M) in ( |20[ ) 
depends on power-spectrum. Since, in general, disjoint vol- 
umes are not independent, this means, of course, that there 
is no longer any guarantee that the algorithm will generate 
the excursion set mean values, either when used as a one 
step partition algorithm, or when the partition algorithm 
has been embedded in a loop over timesteps. Another way 
to see this is to note that f{m\AI) factors conveniently only 
for a white-noise power-spectrum. Since there is no guaran- 
tee that the algorithm will produce the excursion set state- 
ments in the mean, we also have no real analytic estimate 
for the higher order moments associated with this algorithm. 
Later, when we show the results using either the ensemble 
of one-step partitions or multi-step merger trees associated 
with this approach, we will compare the Monte-Carlo results 
with the formulae given in the previous subsection (^, pl| 
and [isl), but with the relation ct^(M) depending on power 
spectrum. 



3.3 Relation to binary split models 

Binary split algorithms based on the excursion set approach 
have been discussed by Lacey & Cole (1993) and by Sheth & 
Pitman (1997). In the limit 5i — 5q = Sio <C 1, the one-step 
partition algorithm discussed above will produce partitions 
in which most of the mass is in the largest two pieces. To see 
this, recall that each choice in the partition process is of some 
s-S{R) = {5io/xf, where a; is a Gaussian random variable 
with zero mean and unit variance. Since x is Gaussian, most 
of time, s - S{R) ^ Sfo- If Sio < 1, then s « S{R) after 
most draws of the random number x. This means that, most 
of the time, the first choice is of a subhalo with most of the 
mass of M. Iterating this argument, the second choice is of 
a subhalo which contains most of the remaining mass. Since 
a given range in s — 5* correponds to a smaller mass range 
for white-noise than for power-spectra with more negative 
slopes, the above is more true for white-noise than for these 
other power-spectra. Therefore, at least for white-noise, a 
binary split algorithm, with f{m\AI) for the split rule, and 
with succesive splits spaced at AS = Sio ^ 1, should provide 
a good approximation to the merger trees produced here. A 
different approach to just such a binary split algorithm is 
described by Sheth & Pitman (1997). They also discuss why 
binary split algorithms that are consistent with the excur- 
sion set approach are easier to construct in the white-noise 
case than in the general case. 



4 MONTE-CARLO REALIZATIONS OF 
MERGER HISTORY TREES 

In this section, the two partition algorithms discussed in the 
previous section are used to produce ensembles of partitions. 
These partitions are compared with the analytic results de- 
scribed earlier. The second algorithm is then embedded into 
a loop over small timesteps to generate ensembles of merger 
history trees. Various properties of these trees are also com- 
pared with the analytic formulae. Comparison with halos in 
N-body simulations is the subject of the next section. 

We show results for partitions associated with initially 
scale free power-spectra: P{k) oc fc". In this case, a'^{m) oc 
m.~", where a = {n + 3)/3. For such power-spectra, ex- 
pressions like N{'m\M) and c{mm\M) are functions of two 
parameters only: m/M and {Si — So) /a{M). It is conve- 
nient to define S = 5c(l + z) = a(A/»)(l + z), so that 
{Si - So)/o{M) = {zi - zo) (M/M.)"''^ The curves below 
are plotted as functions of {m/M)°'. They were produced 
for a range of choices for {zi — zo), with M set equal to M, . 
The scaling above shows that by rescaling {zi — zo) appro- 
priately, the same curves represent other values of M/Mt,. 
For example, if M = fiA'h , then these same curves represent 
the subclump distribution at {zi — zo)/m"''^- 

Figs. |l| and ^ show results for initial power spectra with 
slopes n = 0, — 1, and —2. In all panels, the thick solid curves 
show the analytic formulae with Ad = M, and {zi — zo) = 
0.01, 0.03, 0.1, 0.3, 1 and 3. Symbols show the correspond- 
ing quantities, averaged over 4,000 random partitions of M , 
generated using the algorithms described above. The dif- 
ferent panels in each figure show the cumulative quantities 
F{> m\M), N{> m|M), and Var(> m\A4)/N{> m\M). In 
both figures, F{> m\M) and A'^(> mjM) are got by inte- 
grating the excursion set equations (fed) and (|2l|) from m to 
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Figure 1. Results from using the first algorithm described in the text to partition an M = Af, halo into subclumps at an earlier time, 
specified by (zi — zq) = 0.01. 0.03, 0.1, 0.3, 1 and 3. In the top and middle panels, curves that are higher on the right correspond to 
smaller values of ^51 — zn}. The trend is the opposite in the bottom panel. Symbols show the Monte-Carlo results, curves show analytic 
results (from|2o[ gjand g^. 



M. Var(> m|Af )/A''(> m| A'f) is the variance in the number 
of subclumps more massive than m, divided by the mean, 
A'^(> m\M); if the number of subclumps were Poisson dis- 
tributed among the halos, this quantity would be unity. The 
solid curve for Var(> m\M)/N(> m\M) in Fig. |^ is got by 
integrating equation (M) over the relevant range in mi and 
m2, doing the necessary integral to get the mean A''(>mjA'/), 
so computing the variance, and then dividing by the mean. 



The corresponding curve in Fig. g is got by integrating equa- 
tion (instead of over the necessary range in m. 

The symbols in Fig. |l| are well fit by our analytic formu- 
lae. This shows that the algorithm works as planned. The 
fits to the mean excursion set quantities (top two panels) are 
excellent, as is the fit to the second order moment (bottom 
panel). 

The ratio of the variance to the mean almost always 
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Figure 2. Results from using the second algorithm described in the text to partition an M = M* halo into subclumps at an earlier time, 
specified by (21 — 20) = 0.01. 0.03, 0.1, 0.3, 1 and 3. In the top and middle panels, curves that are higher on the right correspond to 
smaller values of (21 — 20). The trend is the opposite in the bottom panel. Symbols show the Monte-Carlo results, curves show analytic 
formulae (from ^o| ^ and for comparison. 



has a minimum — this is a consequence of mass conserva- 
tion. Consider the range m/M > 0.5 for a given redshift 
interval z = z\ — zq. Let p(n, 2) denote the probability that 
M has n subclumps in this range. Mass conservation requires 
that there be at most one subclump in this range. There- 
fore p{n,z) = for all n > 1, so the mean is /xi = p{l,z), 
the variance is /i2 — p(l,z) — p(l, z)'^ , and the ratio /i2//ii 
is 1 —p{l,z), which is always less than unity. Since p(l,z) 



refers to all subhalos that are more massive than m, it will 
certainly not decrease as m/M decreases from 1 to 0.5 (since 
decreasing m/M corresponds to increasing the allowed range 
of subhalo masses) . Therefore, over this range, /i2 / /ii will de- 
crease as m/M decreases. The rate of decrease will depend 
on the redshift interval, since, for a given m/M > 0.5, p(l, z) 
will decrease as z increases. (This just reflects the fact that 
it is much less likely that m/M « 1 at, say, z = 5 than at, 



© 0000 RAS, MNRAS 000, 000-000 



8 R. K. Sheth & G. Lemson 



^ 0.5 



A 



A 



0.2 

0.1 

10 
1 

0.1 

0.01 

1 



A 0.1 r 



> 



0.01 




10 0.01 0.1 
m/M 



0.1 1 0.2 0.5 1 

m/M) {tu/M) 



Figure 3. Results from embedding the second algorithm described in the text into a loop over small time-steps (A^ = 0.01) to generate 
the merger history tree of an M = M» halo. Figure shows quantities at an earlier time, specified by (21 — 20) = 0.03, 0.1, 0.3, 1, and 3. 
Symbols show the Monte-Carlo results, curves were made similarly to those in Fig. ^ 



say, z = 0.5.) When m/M <^ 1, then the restriction implied 
by mass conservation becomes less severe, so the scatter in 
the subhalo counts can increase, and may even exceed the 
Poisson value. Fig. ^ shows that, for white noise, this Pois- 
son value is never exceeded. The ratio sometimes exceeds 
unity for other power-spectra, at small masses. This is be- 
cause, for this algorithm, subclumps of mass jj, always come 
in groups of > 1, so it is as though less massive clumps are 



clustered, and so the variance in the counts of these clumps 
exceeds the mean. 

The symbols in Fig. ^ show the corresponding results 
generated using our second algorithm (note that the range 
shown on the y-axis in the bottom panel is not the same as 
in the previous figure). For n = 0, of course, there is no dif- 
ference between the two algorithms. For the other values of 
n, the algorithm does not reproduce the excursion set mean 
values exactly, though for m/M > 0.01 or so the agree- 
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Figure 4. The scaled distribution of formation times. Open circles show the distribution generated using our merger tree algorithm, 
solid lines show equation (^), originally derived by Lacey & Cole (1993). The solid lines provide good fits to the distribution measured 
in numerical simulations (Lacey & Cole 1994). 



ment bewteen the symbols and the curves is fair. Over this 
mass range, the variance given by equation ( p3| ) provides a 
fairly good description of the variance in the Monte-Carlo 
ensemble. The algorithm generates partitions with slightly 
too many massive subclumps relative to less massive sub- 
clumps. For n 7^ 0, this means that the discrepancy between 
the algorithm and the theoretical curves is more obvious in 
the plots of F{>m\M) than in those of N{>m\M). 

Fig. ^ shows the result of generating an ensemble of re- 
alizations of the merger history tree of an M = M, halo 
by embedding the second of our algorithms into a loop over 
small time steps {Az — 0.01). The figure shows the same cu- 
mulative distributions as before. Symbols show the Monte- 
Carlo results, and curves were made similarly to those in 
Fig-i 

For n — 0, the symbols are well fit by the curves, and 
the quality of the fit does not depend on Az. This has the 
following consequence. Recall that the first of our algorithms 
uses the white-noise partition of M into subregions, and then 
scales the number of subclumps associated with each subre- 
gion by some factor that depends on power spectrum. Since 
the white-noise algorithm works even when embedded in a 
timeloop, with appropriate care, our first algorithm, when 
embedded in a loop over time steps, will continue to be de- 
scribed by the same formulae that were used in Fig. |l| There- 
fore, we have not included a figure showing this explicitly. 

Fig. ^ shows that when n 7^ 0, merger trees gener- 
ated using the second algorithm are not well fit by the 
analytic formulae. This is hardly surprising, since the one- 
step partitions did not fit these formulae either, and the 
tree is made by stringing together many one-step partitions. 
Since the difference between the symbols and the curves in 
Fig. ^ depended on {zi — zq), we might reasonably expect 



that the Monte-Carlo symbols in Fig. g will depend on the 
choice of Az. We have explored this possibility: for the range 
0.001 < Az < 0.05, this dependence is negligible. This is re- 
assuring, particularly because this means that, at least over 
the mass range m/M > 0.01 or so, equation (^3|) provides 
a good analytic approximation to the true variance of the 
subclump distribution. 



5 COMPARISON WITH NUMERICAL 
SIMULATIONS 

This section compares the ensemble of trees generated by 
embedding the second of our two partition algorithms in a 
loop over time steps with that measured in numerical simu- 
lations of clustering from scale-free initial conditions. First 
we argue that the distribution of halo formation times gen- 
erated by our merger tree algorithm is in reasonable agree- 
ment with what was measured by Lacey & Cole (1994) in 
their simulations. Then we consider a variety of other sta- 
tistical tests, and compare the results generated using our 
merger tree algorithm with those which occured in scale-free 
simulations kindly made available by S. White. 

Suppose we consider the mass of the largest subclump 
of the largest subclump of the largest subclump and so on. 
Define the formation time zt of an Mo halo as the first 
time that this mass drops below half of Mq. Lacey & Cole 
(1993) showed that the distribution of formation times de- 
fined in this way, suitably rescaled, depends only weakly on 
the power spectrum. They showed that, for white-noise ini- 
tial conditions, this scaled distribution is 

dp/dujf = 2 iUf erfc(ajf/\/2), (25) 
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where uj{ = Z( (Mo/M*)"''VV2" - 1. Lacey & Cole (1994) 
showed that this distribution provides a good fit to the 
scaled formation time distribution measured in numerical 
simulations of hierarchical clustering. The solid curves in 
each panel of Fig. ^ show this analytic formula. The open 
circles in Fig. ^ show the corresponding scaled distribution 
of formation times (averaged over 10* realizations of the 
merger tree) generated by our merger tree algorithm. The 
symbols are reasonably well described by the formula. Since 
the formula describes the simulation results reasonably well 
also, we conclude that our merger tree algorithm generates 
formation time distributions that are similar to those in sim- 
ulations. 

Next, we will consider a variety of other statistical tests; 
some are variants of those used by Kauffmann & White 
(1993), and others are variants of those used by Sheth (1996) 
and Sheth & Pitman (1997). The numerical simulations 
used here are the same as those studied by Mo & White 
(1996), where they are described in more detail. They are 
normalized so that the number of particles in an M* halo 
is M«(o) = {a/Sco)^^^"'^^\ where a is the expansion factor 
since the initial time, and n is the slope of the initial power 
spectrum. We show results for n — and n — —1.5 below. 
We will sometimes quote results in terms of M»: for all fig- 
ures below, we use Sco = 1.7, and we study the merger trees 
of Mo halos identified at a time when an M« halo contains 
471 and 166 particles for n = and —1.5, respectively. 

In what follows we choose a minimum mass threshold 
TO. At a given output time we identify all halos with 
mass M > TO. In the remainder of this section, = 0. At 
an earlier output time z\ > we identify all subclumps of 
each halo identified at zo that, at z\ are more massive than 
TO. For each halo, we store the number of such subclumps as 
well as its square. Averaging over all zo halos that have the 
same M/m allows us to compute the mean and the variance 
of the cumulative subclump distribution. In the previous 
section we showed these quantities as functions of m/M as 
TO varied at fixed M, for representative values of M and 
z\ — zo- In Figs. 1^ and ^ we plot these quantities versus 
M/m for fixed m as A/ varies. In all the plots below, when 
zq — 0, then m/M, — 10/471 for the n = simulations, and 
m/M* = 5/167 for n = -1.5. 

Figs. I and | show the first two moments of the sub- 
clump distribution as a function of M/m. The histogram 
shows the number of halos in each mass bin that were 
found in the simulations. The thin solid curves in show the 
expected number computed using the unconditional mass 
function n(M) dM oc f{M)dM/M of equation (|l|). There 
are two reasons for showing this histogram. This first is to 
demonstrate that the unconditional mass function is in rea- 
sonable agreement with the simulations. The second is that 
the histogram is intended to give an idea of the mass range 
over which the halo sample in the simulations is approx- 
imately statistically complete (e.g., some of the statistical 
quantities measured in the simulations may not be accurate 
measures of the true quantities in the mass range where the 
number of halos in a given mass bin is small.) 

The thick solid curves show the conditional mass func- 
tion /ii = N{> TOjM) computed by doing the necessary 
integral over equation (^l|). The thick dashed curves show 
/i2//ii^= Var(> to|M)/A''(> m\M) computed using equa- 
tion (P3), and the thick dot -dashed curves show this same 
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Figure 5. The mean fii, and the variance divided by the mean 
Mz/a*! , of the number of subclumps that, at z, have mass greater 
than m = 10 that are within a halo with mass M/m. Curves show 
the analytic results; for n = these curves describe the results 
of our merger tree algorithm exactly. The filled circles show the 
corresponding quantities measured in the n = simulations. The 
histogram shows the unconditional mass function measured in 
the simulations, and the thin solid curve shows the associated 
Press-Schechter mass function. 
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Figure 6. The same as the previous figure, but for m = 5 and n = 
— 1.5. Filled circles show the simulation results, open circles show 
the corresponding quantities generated using our tree algorithm, 
curves show the analytic formulae, since these formulae do not 
describe the tree results exactly when n ^ 0. 



quantity computed using equation (^^. (For n = these 
two curves overlap exactly.) Recall that the curves in Figs. ^ 
and ^ depend on the combination {zi — zof'S,/S. Although 
the curves there were made for M = M* halos, if M/M, is 
different, then the same curves represent a different value 
of z\ — Zq for scale-free initial conditions. So, with a little 
thought, one can convince oneself that the curves there are 
consistent with those shown here. 

The filled circles show the corresponding quantities 
measured in the numerical simulations. Open circles show 
these quantities generated by our merger history algorithm. 
When n — 0, the tree results are fit by the analytic formu- 
lae exactly, so we only show the theoretical curves. When 
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Figure 7. The distribution of the number of subclumps that, at 
2, have mass greater than m/M* = 10/471 that are within a halo 
with mass Mq/M* = 38/471 = 0.08 (left-most), 447/471 = 0.95 
(middle) and 1024/471 = 2.17 (right-most). Histograms show the 
distribution measured using our algorithm, and filled circles show 
the corresponding distribution measured in the n = simulations. 
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Figure 8. The distribution of the number of subclumps that, 
at z, have mass greater than ra/Mt = 5/166 that are within a 
halo with mass A/q/M* = 28/166 = 0.17 (left- most), 339/166 = 
2.04 (middle) and 1780/166 = 10.7 (right-most). Histograms show 
the distribution measured using our algorithm, and filled circles 
show the corresponding distribution measured in the n = —1.5 
simulations. 



n = 0, these curves fit the simulation results reasonably 
well. When n — —1.5, the merger tree algorithm is not fit ex- 
actly by the analytic curves, and the two analytic curves for 
the variance differ from each other. Figure ^ shows that the 
numerical simulation results are fit better by equation ( |23[ ) 
than by (p^, and still better by the results generated by our 
algorithm. 

Figs. 1^ and ^ show the distribution of the number of 
progenitors of an AIq halo that, at z, are more massive than 
m. The histogram shows this quantity generated by averag- 
ing over 800 realizations of the merger history tree using our 
algorithm, and solid points show the corresponding quantity 
measured in the simulations. For both values of n studied, 
three representative values of Mq are shown: these approxi- 
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Figure 9. The relative sizes of the first and second largest sub- 
clumps. A/i /Mq and M2 / Mi , for the same three values of Mq as 
in Fig. M. Thin histograms show the distribution of Mi/Mq mea- 
sured using our algorithm, and bold histograms show the corre- 
sponding distribution measured in the n = simulations. Open 
circles in the top right show the mean of M2/M1 given Mi/Mq 
generated by our algorithm; filled circles show the corresponding 
mean and rms scatter measured in the simulations. 




0.01 0.1 0.01 0.1 

Figure 10. Same as the previous figure, but for n = —1.5, and 
for the same three choices of Mq shown in Fig. ^ 



mately bracket the range 0.1 < A'l /Af, < 10. For both values 
of n, the histograms generated by our algorithm are in rea- 
sonable agreement with the numerical simulation results. 

Figs.]^ and |l^ show the relative sizes of the subclumps 
produced by our algorithm and those measured in the sim- 
ulations. The histograms show the distribution of the size 
of the largest subclump Mi /Mo for the same values of Mq 
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as in the previous figures. (The topmost panels show the 
lowest values of Mq.) The symbols show the corresponding 
simulation results. The open circles in the top right show 
the mean of M2/M1 given Mi /Mo, where M2 denotes the 
size of the second largest subclump, averaged over 800 real- 
izations of the merger tree. The solid circles and error bars 
show the associated quantities measured in the simulations. 
For both values of n our algorithm is in good agreement 
with the simulations, though the agreement when n = is 
somewhat better than when n = —1.5. 

We conclude this section with the observation that, for 
most of the tests we considered, our algorithm produces 
forests of merger history trees that are reasonably similar 
to those which occur in numerical simulations of hierarchi- 
cal clustering. 



6 DISCUSSION 

We described an algorithm which allowed us to partition ha- 
los at a given time into subhalos at an arbitrary earlier time. 
The algorithm is exact for Poisson and white-noise initial 
conditions. In these cases, we provided analytic expressions 
for the higher order moments of the subclump distribution. 
We discussed two possible ways of modifying the algorithm 
to partition halos which form from more general Gaussian 
initial conditions. We then showed that the one-step par- 
tition algorithm can be embedded into a loop over many 
time steps to generate a forest of merger history trees. This 
also showed why binary split algorithms previously used to 
construct the merger trees may be reasonably accurate. 

We then embedded the second of our partition algo- 
rithms into a loop over time steps to generate the forest 
of merger trees. We compared this forest of trees with the 
ensemble of trees measured in numerical simulations of hier- 
archical gravitational clustering from scale-free initial con- 
ditions. For the initial conditions we studied, the agreement 
between our trees and those of the simulations was fair. Fur- 
thermore, this comparison showed that our analytic formu- 
lae for the higher order moments of the subclump distri- 
bution provided reasonable fits to the numerical simulation 
results even when the initial conditions were different from 
white noise. 

This result is particularly useful for studying the spa- 
tial distribution of dark matter halos. Mo & White (1996) 
argued that the higher order moments of the subclump dis- 
tribution associated with the forest of merger history trees 
of dark matter halos can be related to the higher order mo- 
ments of the spatial distribution of these halos. Since the 
higher order moments associated with our merger tree al- 
gorithm, and also of the merger trees in the numerical sim- 
ulations, were reasonably well fit by our equation (^3|), it 
is possible to write down analytic approximations to spa- 
tial quantities like the halo-halo correlation function that 
are also reasonably accurate. This is the subject of ongoing 
work (Sheth & Lemson 1998). 

Before concluding, we note that Sheth (1998) describes 
a model for clustering from compound Poisson initial condi- 
tions. For these more general initial conditions, the analogue 
of equation (|^) factors similarly to how it does for the Pois- 
son case. This means that a similar algorithm for generating 
the associated forest of merger trees can be used there also. 
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APPENDIX A: THE PARTITION 
PROBABILITY FUNCTION 

This Appendix describes a simple interpretation of the par- 
tition probability function (equation ^). It shows that this 
partition formula is consistent with the assumption that mu- 
tually disconnected volumes in a Poisson distribution are 
mutually independent. 

Since bo can be related to a density (equation ^), an 
(M, 6o)-halo is associated with a region in the initial Poisson 
distribution that has size 

oVM = M/[n{l + So)] = Mbo/n, (Al) 

where n denotes the average density. This means that 
p(ni • • • hm, bi\M, bo) of equation (Q) denotes the probabil- 
ity that a region containing exactly M particles at average 
density n(l-|-5o) has n subregions, each with average density 
n{l + 5i), where 5i > So, of which ni such subregions contain 
only one particle (so they each have size iVi), 712 contain two 
particles (so they each have size 1V2), etc. Notice that the 
sum of the sizes of these n sub- volumes is iVm < oVm , since 
61 < bo- The remaining volume qVm —iVm = M(6o — bi)/n 
is empty. 

Now suppose that the region oVAf containing M parti- 
cles is known to contain n 61-subregions, of which at least 
one such subregion has size -{Vm and contains exactly m par- 
ticles. The probability that oVa/ contains ni 61-regions of size 
iVi, 71,2 61-regions of size 1V2, etc., given that, within oVm, 
there is at least one 61-subregion of size iVm containing m 
particles, and n ~ 1 other 61-subregions is 

(M6oi)"-^e-^^^°^ 77(m,bi) '^rf vii,b,r .^2^ 
v{M,bo) {nm,bi\M,bo} li m\ ' ^ ' 

i — l 

where boi = bo — 61, ^ . = n — 1, and irii = M — m. 
The term 
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{rim, bi nmp{ni ■ ■ ■ nM,bi\M,bo), (A3) 

all parts. 

which is the sum over aU partitions of Ad that have at least 
one m-part, is included in the denominator since this is a 
conditional probability. However, equation (^) shows that 

/ , V ■n(m,bi)'n(M — m,b') 
n„, bi M, bo) = A/601 /,\u, A4 

T](M, bo) 

where b' is defined by the relation 

— =n{l + 5') = -. (A5) 

oVm — iVm 0' 

Thus, b' parameterizes the density in the remaining volume 
oVm — iV m. Substituting equation (A4) for (rim, fcij A/, 60) 
into (1X5) yields 



(A/&01 



\n-2 „-Mboi 



n 



r7(i,&i 



ri{M-m,b') rii! 

i=l 

However, 

{M - m) (b' -bi) = M{bo - fei) = Afboi 



(A6) 



(A7) 



so (A6) can be written as 



[{M - m)ib' - b,)]^-2e-(i'^-^)(^'-'>i) ^^"^ r,ii, bi)"- 

r](M~m,b') 11 ^ 

i = l 

where '^^rii — n— 1, and i = A/ — m. Clearly, this is 
the same as 



p{ni ■ ■ ■ nM-m,bi\M — m, b'), 



(A9) 



where the n^s sum to (n— 1), and the sum of irii is (M — m). 

This shows explicitly that if the region oVm contain- 
ing exactly M particles, so the average density within it is 
M/oVm ~ n{l + So), is known to contain exactly n subre- 
gions, each with average density n(l + 5i), and it is also 
known that one of these subregions has size iVm and con- 
tains exactly m particles, so the remaining M — m particles 
are arranged into n — 1 subregions, each with the same over- 
density Si, in the remaining volume oVm — iVm, then the 
distribution of sizes of the n — 1 subregions (each with over- 
density parameter fei) within oVjvf — iKn is given by the same 
partition probability function as for the total volume (Vm, 
with trivial adjustments for the fact that the average den- 
sity within the remaining volume qVm — iVm will be different 
from that within oVm itself, and for the fact that we need 
now only consider partitions of M—m into n—1 parts, rather 
than of M into n parts. In other words, if p{n,b\\M,bo) is 
viewed as describing the distribution of 61 -subregions within 
the region oVjvf, then it is fully consistent with the fact that 
non-overlapping volumes in a Poisson distribution are mu- 
tually independent. 

This property of the partition formula could have been 
inferred directly from the additive coalescent interpretation 
of equation (^) that is provided by Sheth & Pitman (1997). 



This Appendix shows that the sum over all permutations 
n[m\ of (mi, ■ • ■ , m„), of J^"^]^ 9mi / Ri-i) is unity. This re- 
sult is used in the main text to simplify equation (p^. 
Before proceeding, it is useful to define 



M ■ 



ruk- 



Vi...j...k — ivj — till — ... — iiij 

To illustrate the argument, suppose that n = 3. Since M 
nil + 1712 + ms. 



(BI) 



rrii 



= 1 and 



{rrii + rrij) 



= 1, 



(B2) 



for all combinations of 1 < i,j,k < 3 in which i j ^ 
k. Explicitly, the sum over the n\ — 6 permutations of 
(mi, m2, ms) is 

mi m2 ma mi m2 m2 mi m^ m2 m-j. mi 
M ri ri2 M ri ris M r2 r2i M r2 r23 
^ma mi m2 ^ ms m2 mi 
M rs r3i M rs r32 



The first identity of (B2) shows that this is the same as 



mi m2 mi ma m2 mi m2 ma ma mi ma m2 
M ri M n M r2 M r2 M rs M n ' 
But this is just 

mi (m2 + ma) _^ m2 (mi -I- ma) ^ ma (mi + m2) 



M 



ri 



M 



r2 



M 



and use of the second identity of ( B2 ) reduces this to 



mi m2 ma 

"M "M "M ' 

which is unity, by definition. 

It is easy to see that, for all values of n, the require- 
ment that M = rrii provides relations analogous to those 
in (BS), and that repeated use of these relations will reduce 
the sum of YYi-iij^n/ Ri-i) over all n\ permutations of the 
rrii to unity. 



APPENDIX B: SUMMING OVER 
PERMUTATIONS 

Consider a partition m of M into n parts: m = 
(mi, • • • , m„), where rrii = M. Let Ri = M — X/j=i 
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